module Dqcoeff
using OffsetArrays
using FastGaussQuadrature
export dqcoeff!, dq_coeff, expanded_chebyshev, linspace, logspace, logspace2

function dqcoeff!(xlist::Vector{Float64}, maxorder::Int64, c::OffsetArray{Float64, 3, Array{Float64,3}})
    #xlist: input, one-dimensional real list, with length n
    #maxorder : input, interger
    #c : output, differential quadrature coefficient array, with dimensions (0:maxorder, n, n)
    #implicit none
    #integer(4) :: maxorder, i, j, m, n
    #real(kind=iwp) :: xlist(:)
    #real(kind=iwp) :: c(0:,:,:)
    n = size(xlist,1)
    c .= 0.0
    for i = 1:n
        for j = 1:n
            if i != j
                c[1,i,j] = Mp(xlist,i)/((xlist[i] - xlist[j])*Mp(xlist,j))
            end
        end
    end

    for i = 1:n
        c[1,i,i] = -sum(c[1,i,1:i-1]) - sum(c[1, i, i+1:n])
    end

    for m = 2:maxorder
        for i = 1:n
            for j = 1:n
                if i != j
                    c[m,i,j] = m * ( c[m-1,i,i]*c[1,i,j] - c[m-1,i,j] / ( xlist[i] - xlist[j] ) )
                end
            end
        end
        for i = 1:n
            c[m,i,i] = -sum(c[m,i,1:i-1]) - sum(c[m,i,i+1:n])
        end
    end

    for i=1:n
        c[0,i,i] = 1.0
    end
    return nothing
end

function Mp(xlist, i)
    #implicit none
    #real(kind=iwp) :: Mp, xlist(:)
    #integer(4) :: n, i, j
    n = size(xlist,1)
    r = 1.0
    for j = 1:(i-1)
        r = r * ( xlist[i] - xlist[j] )
    end
    for j = (i+1):n
        r = r * ( xlist[i] - xlist[j] )
    end
    return r
end

function lagrange_basis(xlist0::Vector{Float64}, i::Int64, ξ::Float64)
# xlist0 denotes the expanded-Chebyshev nodes
   len = length(xlist0);
   ξ0 = xlist0[i];
   f = m -> ξ - xlist0[m]
   f1 = m -> ξ0 - xlist0[m]
   iter = [1 : (i-1) ; (i+1) : len]
   u = prod(f, iter)
   v = prod(f1, iter)
   return u/v
end

function dq_coeff(x0::Vector{Float64}, maxorder::Int64=3)
   # x0: expanded Chebyshev nodes, with length n
   # x1: gauss duadrature points, with length n
   ngp = length(x0)
   gx = gausslegendre(ngp)
   x1 = gx[1]
   c = OffsetArray(zeros(Float64, maxorder+1, ngp, ngp), 0:maxorder, 1:ngp, 1:ngp )
   cout = OffsetArray(zeros(Float64, maxorder+1, ngp, ngp), 0:maxorder, 1:ngp, 1:ngp )
   dqcoeff!(x0, maxorder, c)
   for m=0:maxorder
       for p=1:ngp
           for q=1:ngp
               s = sum(c[m, i, q]*lagrange_basis(x0, i, x1[p]) for i=1:ngp)
               cout[m,p,q] = s
           end
       end
   end
   return cout
end

function expanded_chebyshev(ngp::Int64)
    return [cos((2*(ngp - k + 1) - 1.0)*pi/(2*ngp))/cos(pi/(2*ngp)) for k=1:ngp]
end

function linspace(start::Float64,stop::Float64,num::Int64 = 100)
    step = (stop - start)/(num - 1)
    return start : step : stop
end

function logspace(start::Float64,stop::Float64,num::Int64 = 100)
    return 10.0 .^ linspace(start,stop,num)
end

function logspace2(start::Float64,stop::Float64,num::Int64 = 100)
    return 2.0 .^ linspace(start, stop, num)
end

end
